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Abstract. A numerical method for the solution of the elliptic Monge- 
Ampere Partial Differential Equation, with boundary conditions corre- 
sponding to the Optimal Transportation (OT) problem is presented. A 
local representation of the OT boundary conditions is combined with 
a finite difference scheme for the Monge-Ampere equation. Newton's 
method is implemented leading to a fast solver, comparable to solving 
the Laplace equation on the same grid several times. Theoretical justifi- 
cation for the method is given by a convergence proof in the companion 
paper [BF012]. In this paper, the algorithm is modified to a simpler 
compact stencil implementation and details of the implementation are 
given. Solutions are computed with densities supported on non-convex 
and disconnected domains. Computational examples demonstrate ro- 
bust performance on singular solutions and fast computational times. 
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1. Introduction 

The Optimal Transportation (OT) problem is a simply posed mathemat- 
ical problem which dates back more than two centuries. It has recently 
led to significant results in probability, analysis, and Partial Differential 
Equations (PDEs), among other areas. The core theory is by now well es- 
tablished: good presentations are available in the survey [Eva99] and the 
textbook [Vil03]. The subject continues to find new relevance to mathemat- 
ical theory and to applications. 

However, numerical solution techniques for the OT problem remain un- 
derdeveloped, relative to the theory and applications. In this article we 
introduce a numerical method for the OT problem, building on previous 
work, [FOlla, F012] and in conjunction with [BF012]. The solution is 
obtained by solving the Monge-Ampere equation, a fully nonlinear elliptic 
partial differential equation (PDE), with non-standard boundary conditions. 

In [BF012] the convergence proof is established in the setting of viscosity 
solutions. The numerical method is somewhat complicated to explain, so 
this article focuses on implementation details. In addition, we modify the 
algorithm to allow for the use of a compact stencil. While in theory, a wide 
stencil must be used, in our case the compact stencil suffices, except in a 
very singular case which is outside the scope of the convergence proof. This 
allows for a significant simplification of the method. We implement a Newton 
solver, which allows the solution to be computed quickly. Computational 
results are presented including some which are not covered by the theory 
(e.g. non-convex domains, singular Alexandroff solutions). 
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1.1. Optimal Transportation and the Monge- Ampere PDE. The 

OT problem is described as follows. Suppose we are given two probabil- 
ity densities 

px , a probability density supported on X 
pY, a probability density supported on Y 

where X, Y C W 1 are compact. Consider the set, M, of maps which 
rearrange the measure py into the measure px, 

(1) M = {T : X t-> Y, p Y {T) det(VT) = p x }- 
The OT problem, in the case of quadratic costs, is given by 

(2) M±J \\x-T(x)\\ 2 p x (x)dx. 

See Figure 3 for an illustration of the optimal map between ellipses. 
The OT problem (1,2) is well-posed [Bre91]. 

Write Vu for the gradient and D 2 u for the Hessian of the function u. 
The unique minimizing map, M, at which the minimum is reached is the 
gradient of a convex function u : X C M. d — > R, 

M = Vu, u convex ilct^l, 

which is therefore also unique up to a constant. Formally substituting T = 
Vu into (1) results in the Monge- Ampere PDE 

(MA) det( J D 2 n(x)) = P ^f] , for x G X, 

along with the restriction 

(C) u is convex. 

The PDE (MA) lacks standard boundary conditions. However, it is con- 
strained by the fact that the gradient map takes XtoY, 

(BV2) Vu(X) = Y. 

The condition (BV2) is referred to as the second boundary value problem for 
the Monge- Ampere equation in the literature (see [Urb97]). We sometimes 
use the term OT boundary conditions. 

The numerical approximation of the combined problem (MA, BV2, C) is 
the subject of this work. 

1.2. Contributions of this work. The condition (BV2) is not amenable 
to computation. Instead, we replace (BV2) by a Hamilton-Jacobi equation 
on the boundary: 

(HJ) H(X7u(x)) = dist(V«(a;), dY) = 0, for x G dX 

where we use the signed distance function to the boundary of the set Y as 
the Hamiltonian H. 

The problem thus becomes (MA, HJ, C), and solving this system will re- 
quire theoretical and numerical ideas, some of which are addressed in [BF012] 
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The key result from that paper is that H can be discretized using informa- 
tion inside the domain, which results in a scheme that is both consistent and 
monotone on the boundary. Together with a convergent scheme for (MA) 
inside the domain, we proved a convergence result. 

Theorem (Convergence [BF012] ). Let u be the unique convex viscosity 
solution of (MA, HJ, C). Let u e be a solution of the finite difference scheme. 
Then u e converges uniformly to u as e — > 0. 

In the theorem above, the solution depends on a number of parameters, 
which include the grid resolution, the stencil width. These parameters are 
explained below. 

In this work, we describe a numerical method that fits into the conver- 
gence framework of our earlier work. In particular, we describe a method 
for characterizing the target set using a finite collection of boundary points, 
which is the type of information we would typically be given in applica- 
tions. The general theory of [BF012] requires a wide stencil finite difference 
scheme. In this work, we simplify that requirement, and describe a com- 
pact, narrow-stencil version of the scheme, which will nevertheless allow 
us to achieve good accuracy over a wide range of examples from smooth 
to moderately singular. We also discuss solution methods for the scheme, 
and demonstrate how the Monge-Ampere equation and implicit (nonlin- 
ear) boundary conditions can be handled together using Newton's method. 
Finally, we provide extensive computational results that demonstrate the 
capacities of the method. 

1.3. Relation to our previous work. This article builds on a series of 
papers which have developed solution methods for the Monge-Ampere equa- 
tion. The definition of elliptic difference schemes was presented in [Obe06] , 
which laid the foundation for the schemes that followed. 

The first convergent scheme for the Monge-Ampere equation was built 
in [Obe08] ; this was restricted to two dimensions and also to a slow iterative 
solver. Implicit solution methods were first developed in [BFO10], where it 
was demonstrated that the use of non-convergent schemes led to slow solvers 
for singular solutions. In [FOlla] a new discretization was presented, which 
generalized to three and higher dimensions; this also led to a fast Newton 
solver. The convergent discretizations used a wide stencil scheme, which led 
to accuracy less than second order in space. While this cannot be avoided 
on singular solutions, it is desirable to have a more accurate solver on (rare) 
smooth solutions. 

In [FOllb] a hybrid solver was built, which combined the advantages of 
accuracy in smooth regions, and robustness (convergence and stability) near 
singularities. However, this was accomplished at the expense of a conver- 
gence proof. In addition, it required a priori knowledge of singularities, 
which is not available in the Optimal Transportation setting. 

In [Frol2], a more general heuristic method is proposed, consisting in 
iteratively solving (MA) with Neumann boundary conditions, and projecting 
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the resulting set onto the target set Y. The new projection is then used to 
derive new Neumann boundary conditions. This method required several 
iterations, and no convergence proof was available. However, as we note 
in subsection 4.2, this approach can be viewed as a method of solving the 
discretized equations we will describe in this paper. 

In [F012] the convergence theory of Barles and Souganidis was extended 
to nearly monotone schemes, which provided a convergence proof for filtered 
schemes. These approximations retain the stability of the monotone scheme, 
but allow for greater accuracy. In the smooth case, they also remove the 
requirement for wide stencils to be used. 

1.4. Applications. The Optimal Transportation problem has applications 

to image registration [HZTA04], mesh generation [BW09], reflector design [GO03a], 
astrophysics (estimating the shape of the early universe) [FMMS02], and 
meteorology [CNP91], among others. See the recent textbook [Vil03] for a 
discussion of the theory and a bibliography. 

The OT problem also has connections with other areas of mathematics. 
A large class of nonlinear continuity equations with confinement and/or 
possibly non local interaction potentials can be considered as semi-discrete 
gradient flows, known as JKO gradient flows [JK098, OttOl], with respect 
to the Euclidean Wasserstein distance. The distance is the value function 
of the Optimal transportation problem. The impediment so for has been 
the cost of numerical implementation. In one dimension the problem is 
trivial and [KW99] implements JKO gradient flow simulations for nonlinear 
diffusion. An interesting recent work [CM10] considers the two dimensional 
case. The performance of our solver offers opportunities for implementing 
JKO gradient flows. 

1.5. Discussion of numerical methods for Optimal Transportation. 

The numerical solution of the optimal transportation problem remains a 
challenging problem. The early work of Pogorelov [Pog94] introduced a 
constructive method for solving the problem when the target density py 
is in the form of a (possibly weighted) sum of Dirac masses. The convex 
potentials are constructed as a supremum of affine functions whose gradients 
take values in the finite set of supporting points. This method was used to 
build an early numerical solution in [OP88] for a very small sized problem. 
An algorithm for computing these solutions based on lifting hyperplanes and 
Voronoi diagrams can be found in [GM96]. A similar idea has been pursued 
numerically in the contexts of meteorology [CP84], antenna design [GO03b], 
and more recently by in image processing [Merll]. In this context, optimal 
mass transfer is a linear programming problem. 

When the initial density is also a sum of Diracs, the popular auction 
algorithm proposed by Bertsekas (see the survey paper [Ber92]) solves it with 
0(N 2 log N) complexity. In [BoslO], the author compares different linear 
programming approaches and discusses the non-trivial issue of quantization 
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(discretization of densities towards sum of Diracs), which is necessary to 
treat more general mass transfer problems. 

When the density is not simply a sum of a few Dirac masses, other 
techniques are needed for devising computationally practical solution meth- 
ods. One convergent method is the Computational Fluid Dynamics ap- 
proach [BBOO], which relaxes the nonlinearity of the constraints at the cost 
of an additional virtual time dimension. A more recent implementation 
is [HRT10] . 

In two dimensions, complex variables allows the Monge- Ampere equation 
to be solved using a combined numeral analytical techniques, which requires 
the solution of nonlinear equations on the boundary. This technique was 
used to build coarse numerical solutions in [KS84]. 

Most other previous works on numerical optimal mass transportation have 
restricted their attention to simple geometries where a simpler boundary 
condition can be used. One simplification is to work on the torus (with 
periodic densities) using the change of variable u = Id + v (v is peri- 
odic) [LR05, SAK10]. However, this severely restricts the range of solvable 
problems since mass transfer between periodic cells may be optimal. There 
is no easy way to prevent it and capture the optimal transportation for the 
one cell non-periodic problem. 

A second option for designing optimal transportation boundary conditions 
is to consider only very simple geometries, for example mapping a square to 
a square [HZTA04, DCF + 08]. In this case, a Neumann boundary condition 
can be used to produce a face to face mapping on the boundary, which is 
optimal when mass does not vanish. 

1.6. Rigorous approximations of the OT problem. While there are 
a number of implementations of numerical methods for the OT problem, 
many do not focus on mathematical justification for the method. Any prov- 
ably convergent approach the solution of the OT problem in the continuous 
setting must use an appropriate notion of solution. To build a convergent 
method in a general setting , some notion of weak solutions must be used: 
previous experience with the simpler case of Dirichlet boundary conditions 
shows that even mild singularities can affect solver performance dramati- 
cally [BFO10]. 

Solutions of the OT problem lead to singular solutions of (MA, HJ, C), so 
some notion of weak solutions is needed. Various notions of weak solutions 
are available: Pogorelov solutions, Alexsandrov solutions, Brenier solutions, 
and viscosity solutions, see [Vil03, Chapter 4] for references. The first three 
are specific to the OT problem; viscosity solutions are defined for a wide 
class of elliptic PDEs [CIL92], these are also the least general. 

Remark. We choose to use viscosity solutions, because they are well-suited 
to finite difference methods and because there is a robust convergence the- 
ory [BS91]. In addition, in previous work we have build numerical methods 
for the Monge-Ampere equation with standard boundary conditions (see 
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subsection 1.3). Rigorous implementations using other notions of solution 
are also possible, but so far has not led to an effective numerical method. 



We now describe the numerical approach to solving the second boundary 
value problem (BV2). We implement the boundary condition using the 
signed distance function to the boundary of the set Y. We are able to treat 
the boundary condition using a monotone finite difference scheme, which is 
consistent with the treatment of the PDE (MA). 

The resulting boundary conditions are implicit, in contrast to, for exam- 
ple, Dirichlet or Neumann boundary conditions, which have been previously 
implemented for (MA). However, at each grid point the boundary condi- 
tion can be treated as an implicit equation in a similar manner to how the 
interior grid points are treated: in our case, with a Newton solver. 

2.1. Representation of the distance function. In our implementation, 
the function that defines the boundary of the target is the signed distance 
function 



Then at boundary points x £ dX, the solution to the second boundary value 
problem will satisfy H(Vu(x)) = 0. 

In [BF012], it was shown that this PDE can be written as the supremum 
of linear advection equations. Indeed, this characterization is equivalent to 
representing the convex target set by its supporting hyperplanes. Upwind 
discretizations of this type of equation are standard. However, because the 
equation is posed on the boundary, some of the values required by the dis- 
cretization could lie outside the domain. By relying on the convexity of the 
signed distance function H and the solution u, we simplified the expres- 
sion for the boundary condition, which enables us to describe a convergent 
upwind discretization that relies only on values inside the domain. 

Lemma 2.1 (Representation of distance function [BF012]). Let u 6 C 2 (X) 
be a convex function whose gradient maps the convex set X onto the convex 
set Y : Vu(X) = Y . Let n x be the unit outward normal to X at x 6 dX . 
Then the signed distance function to Y at the boundary point Vu(x) E dY 
can be written as 



2. Approximation of the boundary condition 



(3) 




+ dist(y, dY), y outside of Y 
— dist(y, dY), y inside Y . 





M=i 



where the Legendre-Fenchel transform is explicitly given by 



(5) 




sup {y ■ n}. 

y £dY 
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2.2. Representation of the target set. In the above equation, the geom- 
etry of the convex target set is encoded in the values of H*(n) for each unit 
vector n. In our implementation of the finite difference solver, these values 
will either be given or pre-computed, so there will be no need to re-compute 
these terms each time the signed distance function is evaluated in the solver. 

If the target set Y is represented by its supporting hyperplanes, the val- 
ues of the Legendre-Fenchel transform can be obtained using simple convex 
analysis [BV04]. In practice, however, we expect that the target set will 
be represented by a collection of scattered points on the boundary. In this 
setting, we can approximate the values of H*(n) using (5), where the supre- 
mum will be estimated over the finitely many given boundary points yo- We 
stress again that this only needs to be done once, and the overall computa- 
tional cost of this step is insignificant compared to the cost of solving the 
Monge- Ampere equation as a whole. 

In the implementation, we will further discretize this Hamilton-Jacobi 
equation (4) by computing the supremum over a finite subset of the admis- 
sible directions. In practice, we will simply use a uniform discretization of 
the directions, 



2.3. Discretization of boundary condition (4). In this section we ex- 
plain how to build a correct, monotone discretization of H using points at 
the boundary and on the inside of the domain X. 

For simplicity, we describe the discretization on a square domain. We note 
that a slightly more complicated discretization will generalize this to more 
general triangulated domains. However, by padding the source density px 
with zeros (see subsection 5.1), we can handle different geometries while still 
computing on a simple, square domain. We can also easily generalize the 
discretization to higher dimensions. 

We describe the discretization along the left side of the square domain, 
with normal n x = (—1,0). Along this side, the set of admissible directions 
will be given by 



Then, letting dx denote the spatial resolution of the grid, we can approxi- 
mate the advection terms in (4) by 



This scheme only relies on values inside the square and, because n\ < 0, 
it is monotone. Taking the supremum of these monotone schemes over all 
admissible directions, we preserve monotonicity of the scheme. 

We write a similar scheme on the other sides of the square. At corners, 
we take formal limits of the obliqueness constraint to limit the admissible 



(cos (2nj/N Y ) , sin {2irj/N Y )) , j = 1, . . . 



N Y . 



{n = (ni,ri2) | ni < 0, \\n 



1}- 
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directions to a single quadrant, which ensures that the required information 
will continue to reside inside the square. 



3. Approximation of the Monge- Ampere equation 

Next, we turn our attention to the finite difference approximation of 
the Monge-Ampere equation in the interior of the domain. Because the 
solution of (MA) is unique only up to constants, it is convenient to fix 
a given solution. This can be done using the equation det(D 2 u(x)) = 
Px(x)/py(Vu(x)) + (u). However computing the mean is costly in practice, 
so we replace the mean (u) with the value u(xq) at a fixed grid point. 



3.1. Convergent filtered discretization of the Monge-Ampere equa- 
tion. Our approximation of the Monge-Ampere equation is based on the 
convergent, filtered schemes described in [F012]. This almost monotone 
scheme is related to the work of Abgrall [Abg09] which is a "blending" of a 
monotone and a higher-order scheme. 
These schemes have the form 



In the expression above, MA e M is a monotone scheme, which is convergent 
but typically low-accuracy, and MA e A is a more accurate approximation. 
The perturbation e(h, d9) should go to zero as the grid is refined (that is, as 
the discretization parameters h, d6 — > 0). Here the filter function S is given 



det(D 2 u(x)) 



Px{x) 



+ u(xq) for x £ X. 



p Y {Vu{x)) 




by 



x \\x\\ < 1 



(6) 



S(x) 



= < 



||x|| > 2 

-x + 2 l<x<2 
-x-2 -2<x<-l. 



See Figure 1. 
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y 



1 -- 



y = S(x) 



3 




2 



3 



X 



Figure 1 . Filter function 



In order to build a consistent monotone scheme in general wide stencil 
schemes must be used. However, it is observed in [F012] that for smooth and 
even moderately singular examples, the discretization error of the filtered 
scheme does not depend on the stencil width. Thus, except in the most 
singular cases, we can achieve good accuracy using a compact, narrow stencil 
discretization. Here for simplicity and clarity, we describe the compact 
version of the scheme. This method can be extended to a wider stencil if 
higher accuracy is desired on the most singular examples; see [Frol2]. 



3.2. Finite difference operators. Before we describe the discretization 
we use, we define several finite difference operators that we will utilize for 
approximating the first and second partial derivatives using centred differ- 
ences. 




\Dx2 ^] i j 



1 



(Ui,j + 1 - Uij-l) . 



2dx 



We will also be making using of derivatives along the directions 
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The finite difference approximations of these derivatives are given by 

[V vv u]ij = ^2 + 1H-1J-1 ~ 2u i,j) 

[T> v x v xu]ij = ^2 + u i+hj-i ~ 2u hj) 

\ V vMij = 2 ^2 dx ( U i+hj-l ~ u i-lj+l) ■ 

3.3. Compact monotone discretization. The monotone scheme for the 
Monge-Ampere equation is based on the following variational characteri- 
zation of the Monge-Ampere operator for convex functions [FOlla, Frol2]. 
The convexity constraint is enforced by additionally replacing the directional 
derivatives with their positive part. In addition, to prevent non-convex so- 
lutions when the right-hand side vanishes, we will also subtract the negative 
parts of these second derivatives, 

{d d 
IT m.ax.{u VjVj , 0} + m.m.{u v Uj , 0} 
j=i 

which is valid when u is convex. These modifications ensure that a non- 
convex function cannot solve our Monge-Ampere equation (with non-negative 
right-hand side). 

For clarity and brevity, we describe in detail the compact (low accuracy) 
form of the monotone discretization in two dimensions. However, this can 
also be generalized to higher dimensions and wider stencils. 

In the compact version of the scheme, the minimum in (7) will be ap- 
proximated using only two possible values. The first uses directions aligning 
with the grid axes. 

AfAi[u] = max {£>2 ia;i u, 5} max {T> X2X2 u, 5} 

- min{2?a. ia!1 u,(5} - mm{V X2X2 u, 5} - px/pr (D Xl u,V X2 u) - u . 

For the second possible value, we rotate the axes to align with the corner 
points in the stencil. The resulting approximation of the equation is now 

M^H = max {V vv u, 5} max {V v ± v ±u, 5}— min {T> vv u, 5}— min {V v ± v ±u, 5} 

- — \^-=(V v u + T> xu), -^=(D v u - V ±u) J - u . 
Py VV2 V2 J 

The monotone approximation is then given by 

MA M [u] = min {MA X [u] , M A 2 [u] } . 

The effect of this approximation is to fix the directional resolution (which is 
referred to at d6 in [BF012]) to a fixed value. While d9 is required to go to 
zero for the proof of convergence, even for the filtered scheme, in practice, for 
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the solutions we computed, this is not needed. This results in a significant 
simplification of the resulting scheme. 

To ensure monotonicity of this scheme, the parameter 5, which is used to 
bound the second derivatives away from zero, should satisfy 5 > Kdx/\/2 
where K is the Lipschitz constant (in the y variable) of the right-hand side 
Px{x)/p Y {y). 

3.4. Accurate discretization. The filtered scheme also uses a (formally) 
more accurate discretization of the equation. We apply a standard centered 
difference scheme, using the operators defined in subsection 3.2. 

MA A [u] = V xlXl uV X2X2 u - (D XlX2 u) 2 - px/py (T> xl u,T> X2 u) - u . 

4. Solution methods 

In the previous sections we focused on a convergent discretization of the 
PDE, which in our case results in a system of nonlinear equations. In this 
section we focus on the methods used for solving the resulting system of 
equations. In the first case, convergence refers to the limiting solution as 
the discretization parameters go to zero (and monotonicity is important). 
The solutions methods are usually iterative, and in this context convergence 
refers to the convergence of the iterative method. There are monotone so- 
lutions methods, but these are usually slow. However monotonicity is not 
necessary for the solution method. 

The simplest, but slowest method, is explicit iteration. The fastest method 
is Newton's method, but it involves regularization and setting up the lin- 
earized equations. We also describe the projection method of [Frol2], which 
we interpret as a particular semi-explicit method for the discretized system 
(MA-HJ), even though the current formulation was not available at the time. 

4.1. Explicit iteration. The simplest approach to solving the nearly- monotone 
system is to simply perform a forward Euler iteration on the parabolic equa- 
tion 

Ut = det(D 2 u) = px/pyC^u), x G X. 

In order to enforce the boundary conditions, we can simultaneously evolve 
the Hamilton-Jacobi equation 

u t = H(Vu), x e dX. 

By allowing this system to evolve to steady state, we can obtain the solution 
of the discrete system. 

While this explicit solution method will converge to the solution of the 
discretized equations, it is subject to a very restrictive nonlinear CFL con- 
dition [Obe06]. Consequently, this approach is very slow. 
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4.2. The projection method. A much faster approach, which can be 
viewed as a method of implementing (BV2), is the projection method de- 
scribed in [Prol2]. We show that it can be viewed as a simple iterative 
relaxation strategy for solving (MA-BV2). 

The idea of this approach is that we wish to solve for the solution u to 
the Monge- Ampere equation in the domain, as well as its gradient p = Vu 
on the boundary of the domain. We use a splitting approach that involves 
alternating between: 

(i) the solution of a Monge- Ampere equation with Neumann boundary 
data obtained from the current estimate of the gradient map p at 
the boundary. 

(ii) the solution of the Hamilton-Jacobi equation for the gradient p at 
the boundary. 

More precisely, step (i) involves updating u k+1 by solving the following 
Monge- Ampere equation with Neumann boundary conditions. 

'det(D 2 u k+1 (x)) = Px (x)/p Y (Vu k + 1 (x)) + u k+1 xeX 
< X7u k+1 • n x = pk ■ n x x £ dX 

u k+1 is convex. 

v 

The gradient map X7u k+1 obtained from solving this PDE may not solve 
the correct Hamilton-Jacobi boundary condition. Thus step (ii) involves 
updating the value of p k+1 to ensure that it does satisfy H(p k+1 ) = 0. To 
ensure that we are making use of the information from the Monge- Ampere 
equation, we define an intermediate value of the map p k+1 / 2 = \7« fc+1 at the 
boundary. Next, we update using a gradient descent approach. 

p k+l = p k+l/2 _ tf(pfc+l/2) V #(pk+l/2). 

This step has a geometric interpretation. If Vu fc+1 is not in the target set, 
then the new value of p k+1 is precisely its projection onto the boundary of 
the target (dY). Consequently, we recover the projection method of [Frol2]. 

4.3. Newton's method. Because we have succeeded in rewriting the bound- 
ary condition as a PDE posed on the boundary, we could treat the result- 
ing nonlinear equations implicitly, just as we do with the Monge- Ampere 
equation in the interior. We now describe the use of Newton's method for 
accomplishing this. 

4.3.1. The Monge-Ampere equation. We begin by reviewing the form of 
a damped Newton's method for the filtered discretization of the Monge- 
Ampere equation in the interior of the computational domain. This involves 
performing an iteration of the form 

u k+i = u k_ a (yMA[u k ])- 1 MA[u k ] 

where the Jacobian is given by 

VMA[u] = (1 - S'[u]) VMA M [u] + S'[u]VMA s [u\. 
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The derivative of the filter function, given by (6), is 
S'(x) - 



1 ||x|| < 1 
-1 1 < ||z|| < 2 
INI > 2. 



To prevent this from taking on negative values, which can lead to poorly 
conditioned or ill-posed linear systems, we approximate the Jacobian by 

VMA[u] = (1 - S'[u]) VMA M [u] + max{S"H, 0}VMi s [4 

The damping parameter a is chosen to ensure that the residual is decreas- 
ing. 

This iteration requires the Jacobians of both the monotone and accurate 
schemes. To simplify the expression, we define F(x,p) = px(x)/py{p)- 
We also use lo to denote the matrix that has entries equal to one in the 
column corresponding to the point Xq. We begin with the monotone scheme, 
recalling that this discretization has the form By Danskin's Theorem [Ber03] , 
we can write the Jacobian of this as 



VMA M [u] 



{ VMA 1 [u] if MA M [u] = MAi [u] 
\VMA2[u] otherwise. 



Thus we need to compute the Jacobians of the two component schemes. We 
give the values of the row corresponding to the i th grid point. 

V Ui MA l [u] = nwx{Dx B(Ba tj i ,0}lx». 1 > - l^^^o V XlXl 



+ 
dF 



m&x{D XlXl Ui,0} t Vx 



V, 



dF 



(x,V Xl Ui,T> X2 Ui)V Xl - ^— (x,T> Xl Ui,V X2 Ui)T> X2 - 1 , 

Opi OP2 



+ 



max{V vv Ui,0} t v 
1 



1 



V 



v., 



1 OF ( 1 1 \ 

The Jacobian of the accurate scheme is (in two dimensions) 

y ui MAs[u] = (D X2X2 Ui)V xlxl + (V Xlxl Ui)V X2X2 + 2(V xlX2 Ui)V XlX2 

dF dF 

- -^—(x,'D Xl Ui,V X2 Ui)V Xl - -—(x,V Xl Ui,V X2 Ui)V X2 - 1 . 

Opi Op2 
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4.3.2. The boundary condition. The Newton update is also used to enforce 
the correct boundary conditions. 

We recall that the discrete form of the boundary condition obtained in 
section 3 can be written in compact form 

H[ui] = max |max{ni, 0}"D~Uj + min{ni, OjPjTuj 

neAi 

+ max{n 2 ,0}£>~iii + min{n 2 , 0}£>~Uj - H*(n)} . 

Again we stress that at each point X{ on the source boundary, we only need 
to check the directions n that allow for upwinding. We also recall that these 
directions, and the corresponding values of H*(n), are given or computed 
before the start of the Newton solve. 

Now at boundary points, the Newton step will take the form 

u k+l = u k - a (vHiu^y 1 H[u k \. 

As with the Monge- Ampere operator itself, we can use Danskin's Theorem 
to compute the Jacobian at these points: 

V Ui H[u] = max{ni,0}P~+min{ni,0}P^ + max{n2,0}P~+min{n3,0}Py 

where again, the value of n used is simply the direction that is active in the 
maximum. 

5. Numerical implementation 

In this section the details of the numerical implementation are presented. 

5.1. Extension of densities. When computing with finite difference meth- 
ods, it is most convenient to work in rectangular domains. However, it is 
often desirable to solve the mass transport problem in more general domains. 
A simple solution, allowed by our solver, is to extend the density function px 
into a square, assigning it the value zero at points outside the set X that we 
are interested in. This will lead to a degenerate Monge-Ampere equation, 
but the filtered scheme is robust enough to handle this problem. As we will 
show in the numerics section, this approach allows to treat non-convex or 
non-connected domains, or even Dirac measures in the source density. 

We emphasize that this trick cannot be used to extend the target den- 
sity py into a square. This is because the convergence and discretization 
error of the scheme are dependent on the Lipschitz constant of px(x)/py (y)- 
Even if the target density is smoothly extended to have a small positive den- 
sity e outside the target set Y, this function will have a very large Lipschitz 
constant, which makes it impractical to obtain computational results with 
any reasonable accuracy. 

However, it is still important to extend the target density so that it is 
defined in all space. This is because, while the given optimal transportation 
problem will require a density py that is defined only in the set Y, we may be 
required to compute this density function at other points during the process 



16 



J.-D. BENAMOU, B. D. FROESE, AND A. M. OBERMAN 



of numerically solving the equation. As we have just noted, we cannot simply 
allow the density to vanish outside the target set. Instead, we must use a 
positive, Lipschitz continuous extension of py into all space. In some cases, 
when py is a given function, there is an obvious way of extending it into all 
space. Lipschitz extensions can always be obtained by using, for example, 
the method of [Obe05] . The resulting function p Y (y) can always be bounded 
away from zero by considering max{/5^(y), po} for some < po < min py(y). 

5.2. Sources of error. 

5.2.1. Discretization of Monge- Ampere operator. There are several sources 
of discretization error in this method. The first source is the discretization 
of the Monge-Ampere equation in the interior of the domain. This dis- 
cretization includes three small parameters that contribute to the error: the 
spatial step size dx coming from the number of grid points Nx along each 
dimension, the angular resolution d9 coming from the width of the stencil, 
and the parameter 5 that is used to bound the eigenvalues of the Hessian 
away from zero in section 3.3. 

The precise accuracy we can expect to achieve will depend on the regu- 
larity of the solution. For sufficiently regular solutions (which includes some 
moderately singular solutions), we expect the filtered scheme to reduce to 
the more accurate scheme, which has a formal accuracy of 0(dx 2 ). For 
singular solutions, we cannot expect high accuracy. With this in mind, we 
simply use the narrow (9 point) stencil version of the scheme, which we have 
described in this paper. 

We previously observed that for convergence of general problems, the 
parameter 5 should be O(dx). However, when the target density py is 
constant, 5 can be arbitrarily small. Thus we set the parameter 5 = dx 2 so 
that it will not have an appreciable effect on the formal discretization error. 
The filtered scheme also requires us to define the maximum size of the non- 
monotone perturbation; in our computations we set this to \fdx + dO, with 
d6 = 7r/4 for the compact scheme. 

5.2.2. Discretization of transport equations. The implementation of the bound- 
ary conditions requires the discretization of a number of transport equations. 
We use simple upwinding to accomplish this, which leads to a discretization 
error that is 0{dx). 

5.2.3. Discretization of target set. If the target set is polyhedral, we can 
determine exactly which vectors n are needed to represent the target, as 
well as the exact value of the Legendre-Fenchel transform of the distance 
function. 

For nonpolyhedral targets, we approximate the target set by a polygon, 
which introduces additional error. It turns out that at no significant com- 
putational cost, this error can be made insignificant. In the discretization 
of the boundary condition (5), we will consider only Ny different directions 
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rij = 2nj /Ny, which represent the normal vectors to a polyhedral approx- 
imation. The boundary condition is the supremum of functions that are 
linear in n, so we expect this approximation to lead to a first order error 

0(1/Ny). 

6. Computational examples 

We now provide a number of challenging computational examples to 
demonstrate the correctness and speed of our method. 

Whenever an exact solution is available, we provide the maximum norm of 
the distance between the mappings obtained from the exact and computed 
solutions (u ex and u comp respectively): 

max \\Vu ex (x) - Vu comp (x)\\ 2 . 

We also provide the number of Newton iterations and total computation 
time for computations done using the most refined target boundary (largest 
value of Ny). We note that there is no appreciable difference in computation 
time as the value of Ny is varied from 8 to 256. 

Remark (Reading mappings from the figures). In the figures which follow, 
the image of the cartesian grid on the source domain is plotted in the target 
domain. The optimal mapping can be interpreted from the figures by noting, 
first, that the centre of masses are mapped to each other. Next, moving along 
a grid line in the source, the corresponding point in the target can be found 
by using monotonicity of the map: the corresponding grid point is in the 
same direction. 

6.1. Mapping a square to a square. We begin by recovering a mapping 
with an exact solution, which involves mapping a square onto a square. To 
set up this example, we define the function 

q(z) = {-^ + 7^ + 7^) cos(8vrz) + ^zsmfrz). 

Now we map the density 

f{x 1 ,x 2 ) = 1 + 4(q"(x 1 )q(x 2 ) + q{x 1 )q"{x 2 )) 

+ I6(q(x 1 )q(x2)q"(x 1 )q"(x 2 ) - q'(xi) 2 q'(x 2 ) 2 ) 

in the square (—0.5, 0.5) x (—0.5, 0.5) onto a uniform density in the same 
square. This transport problem has the exact solution 

u Xl (x 1 ,x 2 ) = x 1 + Aq'(xi)q(x 2 ), u X2 (x 1 ,x 2 ) = x 2 + Aq(xi)q (x 2 ). 

This gradient map is picture in Figure 2. Since the target set is a square, 
there is no need to discretize its boundary. Results are presented in Table 1. 
As anticipated in subsection 5.2, we observe first-order accuracy. 
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(a) (b) 

Figure 2. (subsection 6.1) (a) A cartesian mesh X and (b) 
its image under the gradient map Vit . 

Nx Maximum Error L 2 Error Iterations CPU Time (s) 



32 


0.0220 


0.0127 


5 


0.3 


64 


0.0110 


0.0064 


9 


1.3 


128 


0.0055 


0.0032 


9 


6.6 


256 


0.0028 


0.0016 


11 


41.6 


362 


0.0020 


0.0011 


13 


101.9 



Table 1. Distance between exact and computed gradient 
maps for map from a square to a square. The number of 
Newton iterations and computation time are also given. 



6.2. Mapping an ellipse to an ellipse. Next, we consider the problem of 
mapping an ellipse onto an ellipse. To describe the ellipses, we let M x ,M y 
be symmetric positive definite matrices and let B\ be the unit ball in 
Now we take X = M X B\, Y = M y B<i to be ellipses with constant densities 
f,g'm each ellipse. 

In M 2 , the optimal map can be obtained explicitly [MO04] from 

Vu(x) = M y RgM~ 1 x 

where Rq is the rotation matrix 

_ / cos(0) -sin(0) \ 
Ke ~ \ sin(0) cos(#) J ' 

the angle 6 is given by 

tan(0) = trace(M 2 T 1 M 2 7 1 J)/trace(M :! : 1 M J 7 1 ), 
and the matrix J is equal to 

J = Rtt/2 = ( 1 ) ' 
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We use the example 



M x 



0.8 
0.4 



0.6 0.2 
0.2 0.8 



which is pictured in Figure 3. 

Results are presented in Table 2, which demonstrates first order accuracy 
in both Nx and Ny- 



0.5 



-0.5 




0.5 



-0.5 



I 



-1 -0.5 0.5 1 -1 -0.5 0.5 1 



(a) 



(b) 



Figure 3. (a) An ellipse X and (b) its image under the 
gradient map Vu (§6.2). 









Maximum Error 






Iterations 


Time (s) 


Nx 






N Y 












8 


16 


32 


64 


128 


256 






32 


0.1163 


0.0773 


0.0693 


0.0669 


0.0665 


0.0062 


3 


0.6 


64 


0.1188 


0.0403 


0.0302 


0.0291 


0.0282 


0.0283 


4 


1.0 


128 


0.1214 


0.0302 


0.0201 


0.0174 


0.0168 


0.0168 


4 


4.2 


256 


0.1206 


0.0278 


0.0116 


0.0101 


0.0092 


0.0091 


4 


20.8 


362 


0.1175 


0.0291 


0.0098 


0.0063 


0.0057 


0.0056 


5 


43.7 



Table 2. Distance between exact and computed gradient 
maps for map from an ellipse to an ellipse. The number 
of Newton iterations and computation time is given for the 
largest number of directions (256). 



6.3. Mapping from a disconnected region. We now consider the prob- 
lem of mapping the two half-circles 

X = {(xi,x 2 ) | xi < -0.2, (xi + 0.2) 2 + x\ < 0.85 2 } 

U {(xx,x 2 ) I xi > 0.1, (xi - 0.1) 2 + x\ < 0.85 2 } 

onto the circle 

Y = {(yi,y 2 ) I y\ + yl < 0.85 2 }. 
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This example, which is pictured in Figure 4, is degenerate since the domain 
is not simply connected. Nevertheless, our method correctly computes the 
optimal map, as the results of Table 3 verify. 




0.5 



-0.5 




-1 -0.5 0.5 1 

(b) 



Figure 4. (a) Two half-circles X and (b) its image under 
the gradient map Vu (§6.3). 









Maximum Error 






Iterations 


Time (s) 


N x 






N Y 












8 


16 


32 64 


128 


256 






32 


0.0453 


0.0267 


0.0255 0.0258 


0.0259 


0.0258 


4 


0.2 


64 


0.0397 


0.0184 


0.0158 0.0146 


0.0144 


0.0139 


4 


1,2 


128 


0.0392 


0.0097 


0.0063 0.0066 


0.0065 


0.0064 


5 


4.5 


256 


0.0432 


0.0110 


0.0084 0.0087 


0.0086 


0.0073 


5 


24.9 


362 


0.0448 


0.0130 


0.0070 0.0047 


0.0045 


0.0039 


5 


45.0 



Table 3. Distance between exact and computed gradient 
maps for map from two semi-circles to a circle. The number 
of Newton iterations and computation time is given for the 
largest number of directions (256). 



6.4. Inverse mapping. Next, we show that we can use our method to 
recover inverse mappings. In this particular example, we compute in the 
unit square using variable densities in both the source and the target set. 
The target density is simply a gaussian in the center of the domain: 

1 ( 0.5|x| 2 
p Y (y) = 2 + — ^ exp 



0.2 2 V 0.2 2 

For the source density, we use four gaussians centered at the four corners of 
the domain. For example, in the quadrant [—1,0] x [—1,0], we use 

1 ( 0.5b - (-l.-l)l 2 
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These density functions are pictured in Figure 5. 

To visualize the mapping between these non-constant densities, we plot 
several marker curves (in this case level sets) in the domain X, as well as the 
image of these curves under the gradient map. These are also in Figure 5. 

The optimal mapping is computed in two different ways: 

(1) By solving the problem directly. 

(2) By solving the inverse problem (mapping py to px), then inverting 
the resulting gradient map. 

In order to check that these two approaches produce the same result, we 
look at the distance between the two maps max x( zx ||Vu(Vu*(x)) — x\\ in 
(Table 4). Even for this challenging example, which involves either splitting 
the gaussian into several pieces or joining these pieces back together, the 
difference between the two computed maps depends linearly on the spatial 
resolution h. 




-1 -0.5 0.5 1 -1 -0.5 0.5 



(c) (d) 

FIGURE 5. (a) The source density px and (b) target den- 
sity py- (c) Curves in the domain X and (d) their image 
under the gradient map. 

6.5. Mapping from a non-convex source. As another challenging com- 
putational example, we consider the problem of mapping from a non-convex 
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N x 


Max Distance 


Iterations 


CPU Time (s) 






Forward 


Inverse 


Forward Inverse 


32 


0.2337 


6 


3 


0.3 0.2 


G4 


0.1252 


6 


3 


1.1 1.1 


128 


0.0649 


8 


3 


6.3 5.2 


256 


0.0329 


9 


4 


35.8 32.6 


362 


0.0233 


11 


4 


95.2 52.2 



Table 4. Distance between forward map and the inverse of 
the computed inverse map. The number of Newton iterations 
and computation time is also given for the two different ap- 
proaches for computing the map. 



domain, which can lead to a breakdown in regularity. In this example, we 
choose a domain shaped like the letter "C" , which is mapped into the unit 
circle. See Figure 6 for images of these sets, as well as the computed gradient 
map. 




-1 -0.5 -1 -0.5 0.5 1 

(a) (b) 




-1 ^05 6 -1 -0.5 0.5 



(c) (d) 

Figure 6. The boundaries of the (a) source X and (b) tar- 
get Y sets, (c) Marker curves in the domain and (d) their 
image under the gradient map. 
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6.6. Pogorelov solutions. Finally, we demonstrate the robustness of the 
solver when dealing with very singular solutions. The source density is 
taken to be the Lebesgue measure on a closed convex subset of X, but a 
non-vanishing smooth probability density could be used. The target density 
is a weighted sum of Dirac masses, 

j=N d 

Px = lx, Py = ^2 q i 6 yy 
i=i 

This problem is related to the model used in [FMMS02] to reconstruct the 
early shape of the universe. 

Our method allows for a singular source density, px , but requires smooth 
target densities py ■ Experimentally, the computational time is linear the 
number of grid points, but does not depend on the density px- This means 
that the cost is the same for 1 or 1000 Dirac masses. 

An alternate, but much slower, solution method is to approximate the 
density px to a sum of Diracs. Then the OT problem is reduced to a 
classical assignment problem. This convex optimization problem can be 
solved by the Auction algorithm in 0(N 2 log N) operations, where N is the 
number of points used in the quantization (assuming this is bigger than Nd). 
See [BoslO] for a review and also [Merll] for a multiscale approach. This 
method scales poorly with the number of Dirac masses used. 

Given M, the mass of the target density, we set the weight at grid points 
representing the Dirac masses to jfj^z and elsewhere. We solve the prob- 
lem with the corresponding choice of target and source densities, then re- 
construct the potential u and the cells, {Vj}, from the solution using the 
Legendre transform . 

The Pogorelov solution has the form 

u(x) = max {x ■ yj - vj} . 
j=l,...,N d 

Let u* be the Legendre-Fenchel transform of u, which is also the solution of 
the OT problem with the densities interchanged: mapping py (singular) to 
px (smooth). Let 

Vj = {x £ X, u(x) = x ■ yj — Vj} 

be the support set of the jth affine function in the definition of u(x). Then 
a necessary and sufficient condition that a be a solution to the Optimal 
Transportation problem is that 

(8) \Vj\=qj, j = l,...,N d . 
The optimal {vj} are given by 

(9) Vj = u*( yj ), j = l,...,N d , 



We use the MPT Matlab toolbox http : // control . ee . ethz . ch/~mpt/. 
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Figure 7 shows randomly positioned Dirac masses embedded in a square, 
which are mapped to a uniform density on a ball. We plot the cells, and the 
colormap indicates the error with respect to the optimality condition (8). 
We also show the convex potential, noting that there is a one to one corre- 
spondence between the gradient of the affine facets and the positions of the 
target Dirac masses. We use a 256 x 256 discretization. All computation 
are done in Matlab on a 2.2 GHz intel Core 17 Laptop with 8 GB of RAM. 

Table 5 shows the maximum and L 2 percentage error. We recall that the 
convex potential computed by our method is the dual conjugate of the one 
represented in the figures. This potential will be singular at Dirac locations, 
causing lower accuracy in the values we compute for (9). Note, however, 
that the solution has the correct form of a convex piecewise affine potential. 
We obtain the correct number of cells, and the correct optimal map ordering. 
Thus these solutions could be used for initialization of exact combinatorial 
optimization methods. 

Table 5 also shows that, with this dual approach, the number of Dirac 
masses has no impact on the run time. In Table 6 we look at the run times 
and accuracy for 300 Dirac masses when we increase resolution up to 4 
million points (remember the discretisation of the initial square domain is 
Nx x Nx-) The cost of the method is linear until Nx = 1024 because of 
out of core memory overheads. 

Remark. The computations are performed using a compact filtered scheme. 
While this is enough to observe convergence and satisfactory accuracy for 
smooth or even moderately singular solutions, a wide stencil can become 
necessary for more singular examples. This was seen with the cone solution 
in [F012]. Additionally, these Pogorelov solutions (and the cone solution) 
are not actually viscosity solutions, so they are outside the scope of the 
convergence theory. 

N d L°° Error L 2 err CPU Time (s) 



3 


0.05 


0.02 


10.48 


30 


0.48 


0.21 


10.56 


300 


0.56 


0.18 


11.8 



Table 5. Normalized errors (percentage) for Nx = 256. 



7. Conclusions 

We computed solutions of the Optimal Transportation problem using a 
PDE based algorithm for solving the Monge- Ampere equation with corre- 
sponding OT boundary conditions. Convergence of the method was estab- 
lished in the companion paper [BF012]. The finite difference scheme used 
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Figure 7. The convex potential surface (red) and cell pro- 
jections Vj. The colormap indicates the percentage error in 
cell area, for (a) 3, (b) 30 or (c) 300 randomly positioned 
Dirac masses (stars) mapped to a uniform density on a ball. 

N x L°° Error L 2 err CPU Time (s) 



64 


0.93 


0.28 


1.48 


128 


0.96 


0.24 


2.6 


256 


0.83 


0.21 


11.7 


512 


0.66 


0.20 


46.76 


1024 


0.64 


0.20 


281.53 



Table 6. Run time and error (percentage) for increasing 
grid resolution with 300 Dirac masses. 
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a filtered method, combining an accurate discretization with a monotone, 
lower accuracy wide stencil discretization. A compact form of the monotone 
finite difference discretization produced accurate results for solutions which 
ranged from smooth to moderately singular. 

In order to prove convergence, the target density must be positive and 
Lipschitz continuous and non- vanishing, with support on a convex set. How- 
ever in practise, the implementation is robust to source densities that vanish, 
are singular, or have non-convex support. Since the method is symmetric in 
the choices of target and source densities, this allows one of the densities to 
have these very singular properties, while the other should be have a convex 
support and be bounded away from zero and Lipschitz continuous on the 
support. More singular examples were computed where one density was a 
weighted sum of Dirac masses. 

Both the boundary conditions and the Monge- Ampere equation were be 
treated implicitly using Newton's method. The run time of the algorithm 
is experimentally linear in the number of grid points, independent of the 
particular choice of solution (as long as the assumptions above were met). 

An extension is to generalize the method to three dimensions. A three di- 
mensional implementation of the Monge- Ampere solver for Dirichlet bound- 
ary conditions has already been performed [FOlla]. The extension of the 
Hamilton-Jacobi solver on the boundary is straightforward. 

The linear solver in the Newton iterate could be improved. We used 
the Matlab backslash operator, which is a direct linear solver. An iterative 
method for the linear solver could increase performance. 

The method could be used to compute JKO gradient flows, using the 
linearized equations to estimate the gradient in the gradient flow. 
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